CONFINED ELASTIC CURVES 

PATRICK W. DONDL, LUCA MUGNAI, AND MATTHIAS ROGER 

Abstract. We consider the problem of minimizing Euler's elastica energy for 
f— ^ simple closed curves confined to the unit disk. We approximate a simple closed 

T-H curve by the zero level set of a function with values +1 on the inside and —1 on the 

^D outside of the curve. The outer container now becomes just the domain of the phase 

f^ field. Diffuse approximations of the elastica energy and the curve length are well 

^~> known. Implementing the topological constraint thus becomes the main difficulty 

^ here. We propose a solution based on a diffuse approximation of the winding 

^^ number, present a proof that one can approximate a given sharp interface using 

a sequence of phase fields, and show some numerical results using finite elements 

based on subdivision surfaces. 
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1. Introduction 



Elastic structures confined to a certain volume or area appear in many situations. 
For example inner membranes in biological cells separate an inner region from the rest 
of the cell and consist of an elastic bilayer. The inner structures are confined by the 
outer cell membrane. Since the inner membrane contributes to the biological function 
>• it is advantageous to include a large membrane area in the cell. In two dimensions 

elastic structures confined to a plane ball have been experimentally produced by Boue 
\^ et alii [1] (see also P HE])- They show that with increasing length the structures 

cn become more and more complex. We are considering here the problem corresponding 

•O to a one-dimensional closed elastic wire constrained in a two-dimensional container 

of circular shape. More precisely we consider a wire whose equilibrium (i.e. stress-, 
strain-free) configuration is given by a circle of radius L/2'ir, and we suppose that 
both the friction between the wire and the boundary of the container, and the friction 
between portions of the wire that are in contact are negligible. We are interested in 
finding (stable) shapes of the folded wire constrained in the container. More precisely 
we are interested in those shapes that are obtainable via pure bending deformation 
processes starting from the equilibrium configuration (in particular no stretching is 
allowed). 

We adopt the following mathematical description of the problem. We represent 
the container by the closed unit disk -Bi(O) := {2; G M^ : \z\ < 1} and the folded 
(isotropic) elastic wire by an immersion 7 : 5*]^ — ?■ i?i (0) of the circle S}^ of radius 
L/27r in the unit closed disk. As for the bending elastic energy, we consider the 
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classical Euler's elastica energy associated to the immersion 7. The configurations 
we are interested in correspond to the (local) minimizers of the bending energy among 
the closed curves that are supported in the unit ball, and that can be reached via 
a deformation path that starts from the circle Sj^, and along which the following 
three constraints are fulfilled: the length of the immersed curve remains equal to L 
(so that we exclude stretching of the wire); the elastic (bending) energy is uniformly 
bounded; the immersed curve may have multiple self-contact points, but does not 
have "self-crossings" (as this would correspond to self-interpentration of the wire). 
It turns out (see Section |2]) that the class of immersed curves satisfying the above 
constrains corresponds to the closure (with respect to the W^'^-weak topology) of 
length-preserving diffeomorphisms of S}^ into -Bi(O). In this formulation there are 
several intrinsic difficulties. Minimizers (for large prescribed length) are expected 
to have multiple touching points. Therefore the associated Euler-Lagrange equation 
involves several Lagrange multipliers and an explicit characterization of the class 
of curves in which the minimum is attained is difficult to obtain. Furthermore the 
constraints of being confined to the unit ball and of not developing "self-crossings" 
are difficult to maintain in a steepest descent method. 

In this paper we propose a phase field approximation of the above problem. We 
justify our approach by an asymptotic analysis and investigate the problem by nu- 
merical simulations. As we already remarked above, in the original sharp interface 
formulation admissible configurations correspond to immersions that can be approx- 
imated by a sequence of simple and closed curves. Since simple and closed curves 
bound an inner set we can approximate such sets by smooth fields with values close 
to -|-1 inside and close to —1 outside. Prescribing the confinement condition is now 
rather simple: the outer container just becomes the domain of definition for the ap- 
proximating phase fields and a boundary condition ensures that emerging structures 
do not leave the domain. An approximation of Euler's elastica energy is well known, 
implementing the topological condition thus becomes the main difficulty. One nec- 
cessary condition is that the phase field approximation of the winding number has 
to be close to 27i. We will use a gradient fiow for a relaxed diffuse approximation 
of the elastica energy that includes soft constraints for the prescribed length and for 
the winding number. For "generic situations" we observe that this is sufficient to 
keep the right topology, and avoid that phase interfaces cross transversally. How- 
ever, in general this method does not exclude that a phase disconnects into several 
pieces. To deal with this issue we show that an additional variable can detect multiple 
components and can be used to prevent structures from disconnecting. 

Let us remark that the same phase-fields approximation we use in this paper for the 
bending energy, has been successfully used in similar contexts {e.g. [TOl [11], [121 [3 E]). 
The main differences between our results and the previous literature are on the one 
hand the numerical methods we develop to solve the (diffuse interface) gradient flow, 
on the other hand the inclusion of the topological constraint in the energy. 

The plan of the paper is the following. In Section |2] we discuss the constrained 
minimization problem in its sharp interface formulation. In Section |3] we introduce 



the diffuse interface approximation. In Section |4] we will prove that we can approxi- 
mate a given sharp interface configuration with a sequence of phase fields. In Section 
|5] numerical simulations are presented that show that our approach works reasonably 
well for "well-behaved" initial data. A more exotic example shows that our topolog- 
ical constraint is in general not sufficient to enforce the correct topology for phase 
boundaries. In the Appendix we therefore propose an improved formulation for the 
constrained problem and indicate why this will lead to the correct result. 

2. The sharp interface minimization problem 

We first discuss the minimization problem in its sharp interface formulation. Con- 
sider the unit ball -Bi(O) C M^, a given length constraint L > 0, and define the 
following class of admissable curves 



Mr 



< 7 : [0, L] — 7- -Bi(O), 7 is a closed and simple C^-curve, |7'| = 1 [ ■ (2.1^ 



In particular, elements of M^ can be represented by C^-diffeomorphisms of the stan- 
dard sphere. For 7 G M^ Euler's elastica energy is given by 

S(7) := i'\i\s)\'ds. (2.2) 



^0 
We then consider the constrained minimization problem: find the optimal value 

niL := inf B{^), (2.3) 

and characterize minimal sequences and possible limit points. Since we expect touch- 
ing points for the optimal structures, minimizers will in general not belong to the 
class M^. However, we do obtain the following compactness property. 

Proposition 2.1. Let {^k)kim be a minimal sequence in Ml- Then there exists 
'jeW'^^\[0,L]), such that 

7fc ^ 7 weakly m W'^^'^{\^,L\) (2.4) 

jor a subsequence /c — )■ 00. The curve 7 has the following properties: 

•y is C^ — closed, (2.5) 



7([0,L]) C fii(O), (2.6) 

|y(s)| = 1 for alls e [0,L], (2.7) 

7 can touch the unit circle only tangentially, (2.8) 

7 has no transversal crossings, (2.9) 

where the last property means that 7 can touch itself only tangentially. Furthermore, 
if we denote by Af^ : M^ —)■ {0, 1} the characteristic function of the open subset of 
Bi{0) that is enclosed by 7fc we obtain that 

Xk -^ X strongly m ^^(M^). (2.10) 
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The limit characteristic function X has the following properties: 

X = Xei where E C -Bi(O) is a set of finite perimeter, (2-11) 

d*E C supp(7). (2.12) 

Finally, 7 lies always on the same side of E: after changing the orientation of 'j if 
neccessary, 

ue{x) = ^ 7'(s)^ for n^- almost all X ed*E, (2.13) 

'y{s)=x 

where in the last equation ± denotes the clockwise rotation by tt/2 and ue the inner 
unit normal of E. 

Proof. By the minimizing property of 7^ we have that there exists A > such that 

B{-fk) < A. 

We moreover assume that all 7^ are parametrized by arclength. Since 7 maps to the 
unit ball we therefore have a uniform bound of the sequence (7A:)A:eN in W'^''^{0,L). 



We therefore deduce (2.4) and by Sobolev embedding Theorem that 7 G C^'^/^([0, L]) 
and that 

7fc ^ 7 strongly in C^'''{[0,L]) for all < a < J. (2.14) 



This also implies that (2.5)-(2.9) holds. 



For the inner sets we have a uniform area bound by the confinement constraint 
and a uniform bound on the perimeter, that has length L by the length constraint. 
Therefore X^ is uniformly bounded in BViM."^) und we deduce that for a subsequence 
(2.10) holds and that X is the characteristic set E satisfying (2.11), (2. 12). 



Finally we can orient all 7^ such that 7'(s) equals the inner normal of the set 
that is enclosed by 7^. Then we obtain for any function rj G C^(]R^) from the GauB 



Theorem and by (2.14) that 

UEix) ■r]{x)d\VX\{x) 



X{x)'V ■ rj{x) dx 



lim 

fc— >oo 



lim 

fc— >oo 



Xk{x)V ■ rj{x) dx 



7fc(s)^ ■V{lk{s))ds 



7'(s)^ ■r]{-fk{s))ds, 



which proves (2.13) since rj was arbitrary. 



D 



Proposition 2.1 in particular shows that minimizers of (2.3) belong to the closure 
Ml of Ml with respect to the weak-iy^'^([0, L])-topology. For our purposes, however, 
we only need the following alternative characterization of Ml'- curves in Ml can be 
approximated strongly in W"^'^ by closed simple curves that are strictly contained in 
the unit ball. 



Proposition 2.2. For any 7 G M^ there exists a sequence {'yk)keN of simple closed 
C"^ -curves with 

7fc -)■ 7 as k -^ 00 strongly in iy^'^([0, L]), (2.15) 

|7^|(s) = 1 forallse[0,L], (2.16) 

7fc([0,L]) C i?i(0). (2.17) 

Proof, (i) We first assume that 7([0,L]) C -Bi(O) and therefore that 

6 := dist(7,^i(0)) > 0. 

Repeating the proof of pj Corollary 5.2] under the additional hypothesis that, being 
7 e Ml, 7 is the ly^'^- weak-limit of a sequence of diffeomorphisms of the unit circle 
with equi-bounded "bending-energy" , we obtain the existence of a sequence {ck)ke'N 
of simple closed C^-curves and a sequence {Xk)ken of positive numbers such that, 

Cfc -)■ 7 as fc -)■ 00 strongly in PF^'^([0, L]), (2.18) 

141 (s) = Afc for all s E [0,L]. (2.19) 



From (2.18) it follows that A^ — ?■ 1 as A; — t- 00. Let 

5k := dist(cfc,S'i(0)) 



then (2.18) yields 5^ -^ 5 as k -^ 00. We now define 



1k{s) := T-Cfc(s), s G [0,L] 

and observe that 7^ is a simple closed C^-curve with |7^(s) = 1| for all s G [0,L]. 
Moreover we have 

7fe -^ 7 strongly in W^'\0,L), 

dist(7fc,^i(0)) = 1-i^ ^ 6, 

in particular 7a:([0, L]) C -Bi(O) for k large enough. Therefore (7A;)fceN has all required 
properties. 



(ii) We next consider the general case 7([0,L]) C -Bi(O). First we observe that 
7([0, L]) n -Bi(O) cannot be empty for L > 27r (for L < 2it any minimizing sequence 
converges to a parametrization of the circle with length L). In fact, assume the 
contrary and let (q);^^ be a sequence in Ml approximating 7 weakly in iy^'^([0, L]). 
Then, using GauB Theorem 

L = lim / z/i(s) ■ 7i(s) ds = lim / Xi{x)\/-xdx = 2tt, 

which is a contradiction. Therefore 7 cannot be entirely contained in the unit circle. 
Next consider ^ G C^(-Bi(0)). Then there exists to > and a smooth evolution of 
diffeomorphisms (ft '■ -Bi(O) — )■ i?i(0), t G (—to, to), such that 

d 

gT'^tis) = ^{ft{s)) for all s G [0,L], ipo = Id. 
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Then Cj(s) := ^94(7(3)) defines a smooth evolution of C^-closed W'^'^ curves such 
that all Ct are contained in -Bi(O) and such that q — )■ 7 strongly in iy^'^([0,L]) as 
t — )■ 0. We compute for the length L(t) := f^ \c[{s)\ ds that 

"^ L{t) = - [%"{s)-ai{s))ds (2.20) 

t=o Jo 



dt 



and observe that this expression cannot vanish for all ^ G C^(-Bi(0)) since otherwise 
7([0,L]) n -Bi(O) consists of a collection of straight lines, which contradicts the fact 
that 7 can touch 6*1(0) only tangentially. Therefore we find ^ G C°°(i?i(0)) and 
to > such that the length of q is strictly increasing on [0, to) and such that q — )■ 7 
strongly in 14^^'^([0,L]) as t \ 0. In the following we fix such ^ and to and define 
modified curves with length L, 

lt{s) := JjjM^i^)). se [0,L] 

where a{s) denotes the arclength reparametrization, such that |7^| = 1 holds. Then 
7i is strictly contained in Bi{0). Moreover, we claim that 7^ G M^. In fact let 
(c/)/gN be a sequence in Ml approximating 7 weakly in iy^'^([0, L]). Define curves 
ci,t{^) '■= Vt{ci{s)) according to the variation field ^ fixed above. Then it follows from 



(2.20), the choice of ^, and the weak W"^'"^ convergence of q to 7 that 

dt t=o 
for all / large enough. We then set 



L(/,t) = - / c'i'{s)-aci{s))ds > 
Jo 



as above and obtain that 7;,t G Ml. Moreover we have that •yi^t — ^ 7t as / — )■ 00 
weakly in iy^'^([0, L]), hence ■jt G Ml- Thus, we can apply part (i) and obtain a 
sequence of 'yt,k G Ml that approximates 7t strongly in W"^'^. Taking a diagonal 
sequence proves the claim in the general case. D 

3. The diffuse interface approximation 

Phase field approximations of sharp interface problems are widely used for numer- 
ical simulations and arise from mean field descriptions of phase separation processes 
in various applications. In the following m : Si(0) — )■ M is a smooth function. The 
basis of the phase field formulation is an interfacial energy of the form 

L,{u) ■■= - I (^-\Vu\^ + -W{u))dx. (3.1) 

Here £ > is a small parameter and W denotes the standard quartic double-well 
potential 

W{r) = \{l-rr. (3.2) 



It is well-known [15] that —L^ approximates the curve length functional in the sense 
of Gamma-convergence, where 



Co 



f ^/2W{s)ds. (3.3) 



A phase-field analogue of Euler's elastica energy was already proposed by De Giorgi 
[S]. For the modified version 

B,{u) = - [ ^(-eAu+^W{u)]^ (3.4) 

Co Jb,{o) ^ ^ e J 

the approximation property was proved in two and three dimensions [17] . Moreover, 
following [2] we introduce the diffuse winding number 



T^(u) = — I (-eAu + -W{u))\Vu\. 
Co iBi(o) ^ e / 



(3.5) 



Finally we propose to approximate the constrained minimization problem (2.3) by 
the problem of minimizing 



F,{u) = B,{u)+e-''{L,{u)-L) + e'^ (t,{u) - 27:^ (3.6) 

under the boundary conditions 

u{x) = —1, Vm(x) ■ X = for all |x| = 1, (3.7) 

that prevent diffuse interface from touching the outer container. 



The existence of minimizers for (3.6), (3.7) follows with the direct method of cal 



cuius of variations. Since we are interested in minimizers of the functional B the 
adequate statement regarding the relation between the sharp and diffuse minimiza- 
tion problems would be the Gamma-convergence of J-'s to B. Though we are not able 
to prove such result in full generality, nevertheless we do obtain a compactness result 
and a lower bound estimate in the case of a regular limit point as a consquence of 
dZ] (see also pj). 

Proposition 3.1. Let {u£)s>o be a sequence of smooth functions u^ : -Bi(O) — )■ M that 



satisfy the boundary condition (3.7) and assume that 



sup J-'e{Ue) < CXD. (3.8) 

£>0 

Then there exists a set E C -Bi(O) of finite perimeter such that 

Ue -^ 2Xe - 1 strongly m L\Bi{0)). (3.9) 

Moreover, the diffuse interface measures 

/^^ •= -fllVw.l' + -Wiu,)) dx (3.10) 

Cq \Z £ / 
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converge in measure to a Radon measure /i with support in -Bi(O). If fi is given by a 
curve 7 G M^ in the sense of 

r]{x)dfi{x) = / r]{-f{s))ds for all rj e C^{R^), (3.11) 



then 

Bh) < liminf J;(mJ (3.12) 

holds. 

In general, the limit measure /x will not be given by a curve in Mi but will enjoy 
some weak regularity (being an integral varifold with weak mean curvature in L^). 
We will demonstrate in Appendix |X] that fi can consist of several disjoint curves and 
therfore does not belong to Ml- On the other hand we are mainly interested in the 
numerical simulation of a steepest descent evolution for J-'^ and this in fact works 
sufficiently well. In general a more complex functional is needed, and we propose in 
Appendix \K\ a possible choice. 

4. Construction of recovery sequences 

Whereas we cannot prove that minimizer Ug converge to curves 7 G Ml we can 
show that any such curve 7 can be approximated by a suitable recovery sequence. 
This result also extends to the improved functional we propose in the next section 
and justifies our approximation of the sharp interface minimization problem. We 
first start with the most regular case. 

Lemma 4.1. Let 7 G Ml be given. Then there exists a sequence u^ : -Bi(O) — ?■ [—1, 1] 



of smooth phase fields such that the diffuse interface measures /Xg (as defined in ( 3. 10 ) j 
converge to the measure fi that is 

fie ^ CoWVl (4.1) 

as £ — )■ 0. Furthermore for all e > holds 

L,{u,) = L + R['^\ (4.2) 

%{u,) = 27i + Rf\ (4.3) 

B,{u,) = B{^) + Rf\ (4.4) 

where Re ,Re are exponentially small in e > and Rk is of order 0(5^). In 
particular, 

J'eiUe) ^ i5(7) (4.5) 



for any choice of a, (3 > in (3.6). 



Proof. The construction is standard and uses the optimal profile q for the one- 
dimensional minimisation in the Cahn-Hilliard energy, the signed distance func- 
tion from 7, and an interpolation to the stationary points ±1. To be precise, let 



g : M — 7- (— 1, 1) be the solution of 

-g" + W'{q) = 0, (4.6) 

g(-oo) = -1, g(+oo) = 1, g(0) = 0. (4.7) 

Then 

g'(r) = ^/2Wiqir)) (4.8) 



holds for all r > and with (3.2) we have 

q{r) = tanh(r/y2). (4.9) 

Moreover there exists 6 > such that signed distance function d from 7 (taken 
positive in the region inside of 7) is of class C^. Next fix a smooth symmetric cut-off 
function r] G C°°(]R), 

< rj < 1, ri{r) = 1 for r G [—1, 1], rj^r) = for \r\ > 2, f]' < 0. 

We then define 

2?^ v '2,T' 

Qe{r) ■■= ri{—)q{-)+sgn{t){l-ri{—)) 

and 

u,{x) := qe{d{x)). (4.10) 

Step 1: Consider the parametrization 

V': [0,L) X (-5,5) ^ fii(O), ^(s,t) = 7(s)+tz/(s), (4.11) 

which is injective by the choice of 5 and continuously differentiable with 

detL)^(s,t) = l-te(s). (4.12) 

We then compute that 

Le{Ue) = j j (^-(^,{tf+U¥{q,{t)Yl+t>,{s))dtds 

L rS/2 . . 

q'it/^Y + W{qit/e) ) (1 + tK{s)) dt ds 



J-S/2 ^ ^2 



+ r / (^g;(t)^ + V(g,(t)))(l + t«:(s))dtci.. (4.13) 

Jo JU/2<\t\<S\ V^ £ / 
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By the symmetry of q and (4.8 ), (4.9 ) we obtain for the first integral on the right-hand 
side that 



' q'(t/ef + W{q{t/e) ) (1 + tK{s)) dt ds 



J -6/2 
r6/2e 



e\2 



2L I q{t)^2W{q{t))dt 
Jo 

coL-2L [ ^/2W{r) dr 



«(^) 



CoL-v^Lfl-tanhf^^)") - ^-^ f 1 - tanh^ (^^) 



(4.14) 



Furthermore, using (4.8) again 



2 ,.2t 
6 



1 /2t. .A. 



With this equahty and the symmetry of g^ we calculate for the second integral in 



(4.13) that 



{5/2<\t\<S} 
5 



e 



^-q'^{ty + -W{q,{t))]{l + tK{s))dt 



(4.15) 



Together with (4.14) we obtain (4.2) with 



i?l^) 



-V2l(i -tanhf — —) 
^ ^2^e' 



f.(5 



(\ -tanh^( — - 

+2^/, (i-«y) 2(-f'(T'+7i;"yC+«y)) * 



+ 2L/ (l-,(^))^(l + ,(i)))d, 

2 



which is exponentially small in e > 0. 
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Step 2: Let ^ G C°(5i(0)). We compute that 

f^siO = I I [yAt)' + -W{q,{t)y{^{s)+Hs)){l+tK{s))dtds 

r] [q'{t? + W{q)Y{l{s) + etv{s)){l + etn{s)) dtds 

+ 11 r-q'M + -W{qe{t)\{^{s)+tu{s)){l + tK{s))dtds. 

Jo ^{|<|t|<5} V/: £ / 

(4.16) 
As above we conclude that the second term is expentially small in e > and that 



(4.8) we derive 

L 



Jim/ie(0 = Co / ^i-fis))ds = Cq j ^dH 



n~l 



which proves (4.1). 



Step 3: From (4.8) we obtain, using the shortcuts tj = rj{5 ^2t), q = q{e ^t) etc., that 



-eq': + -W'iq,) 
= (1 - g) (e^V' - ^V(l + g) - %(1 - r;)(l - g)) (4.17) 

is exponentially small in e > 0. For the distance function we have 

{Ad){^{s)+tu{s)) = ^^^ (4.18) 

and for the diffuse mean curvature we obtain 

-eAu, + -W'{us) = -eql + -H^'(g,) + eq^Ad. (4.19) 



Therefore 



/ (-£Am, + V'(m,))|Vm,| (4.20) 

f^ /' ( - ecl'At) + V'(g,(t)) + ^^^Wr^^j'^^Wll + ^'^(^)) ^^^^ 

(^ / K{s)ds\ / eqXtfdt 

+ 2(1 Kis)ds) J^ (^_eq':it) + ^W'{qeit))yq'Mdt. (4.21) 
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Since 7 is closed and simple we have Jg k,{s) ds = 27t. Therefore 

Co 



TeiUe) 



27r 



/ sq'M' dt+— (- eq':it) + -W\qeit))Wit) dt 
Js Co i| V e ) 

(4.22) 



and similarly as above one shows that this term is exponentially small in e > 0. 
Step 4: As above we deduce that 



1 /-^ /-M 



Co Jo J-i ^ 



.g:(t) + V'fe(t))+eg^(t): ""^'^ 



l + tfi;(sl 



'1 + tfi;(s)) dtds 



L rS 



Co Jo Jo 

-L rS 
+ 



^''("^"W'li + M.) 



1 



1 



1 -tK S] 



dtds 



[ [ (-eq'^{t) + -W'{q,{t))y-{l+tK{s))dtds 
Jo Jo ^ ^ ' ^ 







K{s?ds\- f'^ q'{t)^dt 
Co Jo 



Co Jo Jo 1 - eH^K{s)^ 



+ 



Co Jo Jf ^ 
-L rS 



m'<s) 



1 



dtds 
1 



l+tK{s) l-tK[s) 

21 



dtds 



+ 11 (-eq':{t) + ^W'{qe{t))) -{l+tK{s))dtds 
Jo Jo ^ ^ ^ ^ 



(4.23) 



The last two terms on the right-hand side are exponentially small and we finally 
obtain (lO) with Rf^ = 0[e^). 
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We next can prove the general case. 



Proposition 4.2. Let 7 G Mj^. Then the same conclusions as in Lemma 4-i 



hold. 



Proof. By Proposition 2.2 we can approximate 7 strongly in W'^''^{ 0,L]) by a se- 
quence of closed simple C^-curves (7fc)feeN that satisfy (2.16) and (2.17). In particular 
7a: £ Ml and since 7^ — ;■ 7 strongly in W^^'^([0,L]) 

B{jk) ^ i3(7) (4.24) 



as A; — ;■ 00. Lemma 4.1 yields a sequence of functions {us^k)k(^N that satisfy (3.7) and 

J^e{Ue,k) -^ Bi'Jk) 

as e — 7- 0. Choosing now a suitable diagonal sequences proves the claim. D 
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5. Numerical simulations 

In order to demonstrate the feasibility of the above phase field approach to model 
confined elastic curves we present some numerical results. To be exact, we use a finite 
element approach to discretize a viscous gradient flow of the energy T^i after some 
modifications described below, in space and advance the equation in time using a 
first order fully implicit scheme. For other numerical approaches to a diffuse interface 
approximation of constrained Willmore flow see for example p3l ITTl [TO] . 

5.1. Evolution equation in the numerical simulations. As it turns out, for 
finite epsilon, the numerical method does not always yield a perfect transition layer. 
For large prescribed length L it can be energetically favorable to not follow the 
optimal profile of the transition layer — thus increasing the value of the diffuse length 
functional — in cases where two transistion layers were close together. It was therefore 
necessary to introduce a penalty for the discrepancy of the phase field to the optimal 
profile. This term is 



M,{u) = a„i, / (l \Vuf - -W{u)] 



It is evident from the proof of Lemma 4.1 that the addition of such a term does not 



change the construction of the recovery sequence and simply vanishes in the limit of 
small e ii a scales as some power of -. 

Unfortunately, the non-differentiability of the factor |Vu| in the diffuse winding 
number proves to be another problem for the gradient flow. Its gradient yields t^^, 
so the second derivative blows up where |Vm| vanishes. Using the fact that the 
discrepancy of the phase field and the optimal profile have to vanish, we have 



/2 1 

\Vu\ = —^/W{u) = ^ \l 



.21 



U 



The second derivative — which is necessary for the Newton-Raphson iteration used in 
the implicit time integration — of this term still blows up when m = ±1, however, the 
phase field should remain in the interval [0, 1]. For the computation, we thus simply 
leave out the absolute value in this term and observe that the phase field behaves 
nicely in the simulation. 

In conclusion, we numerically compute the viscous gradient flow of the energy 

Te{u) = B,{u) + £-" (L,(n) - Lf + cpe'^ {T,{u) - T)' + M,(n), (5.1) 
where 

Ts{u) = - I (-e^u + -W{u)) ^ (1 - u') . 

The boundary conditions are clamped, i.e., m = — 1 on the boundary of the domain 
and the normal derivative of u on the boundary vanishes. 
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Name 


Length 


Winding Number 


Mismatch 


Length 




Constraint 


Constraint 


Penalty 


Target 


Circle 1 


off 




off 


off 


n/a 


Circle 2 


off 




on 


on 


n/a 


Relaxation 


on 




on 


on 


8.7838 


Topology 1 


on 




off 


on 


8.7838 


Topology 2 


on 




on 


on 


8.7838 



Table 1. The parameters for the numerical experiments 



5.2. Numerical method. The space discretization of (5.1 ) requires some care, since 



its weak formulation requires u to be in H'^{Bi{0)), making it impossible to use a 
piecewise linear interpolation directly. While there are several options to resolve 
this problem, we resort to using a conforming, i.e., continuously differentiable finite 
element discretization. To this end, we construct basis functions derived from Loop 
subdivision surfaces, which can be thought of as a generalization of multivariate 
splines to tessellations of arbitrary topology [13]. The use of subdivision surfaces 
for this problem has been suggested in [7], where one can also find a description of 
convergence properties. In addition, we use the method described in [3] to fix the 
clamped boundary conditions. The computational domain is a disk of radius one, 
discretized using distmesh [16j. In order to advance the system in time, we use a 
simple first order implicit Euler scheme, since accuracy of the time integration is not 
our primary concern. 

5.3. Simulation parameters and results. We use a triangulation of the domain 
consisting of 17813 faces. The transition length e is kept fixed at 0.025, much larger 
values produced a significant mesh effect. The parameters a and (3 are fixed at 2 and 
C/3 = 3. The mismatch penalization o"mis is 0.025~^. 

For the numerical method it is essential to impose initial conditions that already 
are close to an optimal profile of a simple closed curve. To generate such initial 
conditions, we take black and white image to represent the interior and exterior of 
the initial curve, apply a Gaussian blur and use the grayscale data as the initial 
function values. It is then necessary to relax this initial condition in order for it to 
be close enough to an optimal profile for the penalty terms to make sense. We thus, 
in the beginning, chose a small timestep (10~^), and slowly increase the penalization 
of the length- and the topological constraint. The initial conditions are plotted after 
this relaxation phase which lasts 200 timesteps. The timestep is then increased to 
the regular value of 10~^. In addition, we slowly increase the target value L for the 
length constraint, starting at the value of the diffuse length functional at the initial 
condition (after relaxation). For comparison, we also provide some simulation results 
lengh or winding number constraints. 

In the following, we briefly describe the simulation results. The parameters for the 
various simulations are indicated in Table [l} There, "on" for a penalty term means 
that the respective term is used as in equation (5.1). "Off" means the term is not 
present in the energy used for the computation. 
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(a) Time vs. radius plot 




Difference 
0.01 

lo.oos 
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Q 
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-0.01 
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(b) Difference between the expanding cir- 
cle solution with and without winding 
number and mismatch penalization 



Figure 1. Expanding circle 



Expanding circle (Circle 1—2). Figure [I] shows the expansion of a circle of initial 
radius 0.3 with and without length constraint. While it is clear that the phase 
field approximation of Euler's elastica energy B^ alone gives a good approximation 
for the Willmore flow of a radially symmetric initial condition, one can see from 
the radius- vs. -time plot that the topological constraint does not influence this rate of 
expansion. The difference of the phase field of the two simulations (with and without 
topological and mismatch penalty) is shown in Figure 1(b) Note that the maximum 
of the deviation is small compared to 1. 

Relaxation of a folded mirror-symmetric structure (Relaxation). It is clear 
that the gradient flow routine will only find local minima of the energy, and it stands 
to reason that there are many such local minima. We want to investigate the relaxed 



energy of the gradient flow with initial condition shown in Figure 2(a) Figure |2(b) il- 
lustrates the evolving surface. The final relaxed state, with its diffuse energy overlaid, 
can be seen in Figure 2(c)[ The final energy is 33.6. 

Topological transition (Topology 1—2). We investigate the effectiveness of the 
penalization of the diffuse winding number as it differs from 27r. To this end, we 



start a simulation with fairly high energy in the state illustrated in Figure 3 (a 



Figures 3(b) and 3(c) show the state at t = 0.08 for the simulation not penalizing the 



diffuse winding number and penalizing the diffuse winding number, respectively. One 
can clearly see that the simulation without penalization is getting close to pinching 
off at two positions. Finally, one can see that a topological transition occured in 



Figure 3(d), while the curve in Figure 3(e) remained simply connected. Both figures 



are taken at t = 0.1. The overlaid diffuse winding number functional in those figures 
clearly shows how the topological transition changes the calculated winding number 
integral. 
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-OS -0-6 -0 4 -0 2 2 4 06 



(a) Initial condition (after initial 
ation) 



relax- 



(b) Zero level sets of the intermediate 
stages, the arrow indicates the direc- 
tion of the flow. 




Energy 
0.075731 



0.06 
0.04 



-0.02 

t 

6.758e-21 



-08 -06 -0.4 -0.2 0.0 0.2 04 0.6 0.8 



(c) Zero level set of the equilibrium state. 
Elastic energy density is shown as overlay. 

Figure 2. Mirror symmetric initial condition 

Appendix A. A better topological constraint 

The numerical simulations presented in this paper suggest that the "topology" 
of the diffuse interface is preserved along the gradient flow of J-'e when the initial 
condition is well-prepared around an element of Ml and L is not too large with 
respect to the diameter of the domain. However, in general, neither the functional 
J-e nor the winding number in the sharp- interface setting enforce the correct topology, 
as the following example shows. 
Consider 



E = (^5i/2(0) U 5i/5((0, 3/4)) j \ 5i/4(0) (A.l) 

and consider diffuse approximations u^ obtained via the construction presented in 



Lemma 4.1 It is then easy to see that \Te{ue) — 27r| is exponentially small in e. 
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(b) Critical stage of topological transition 
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winding number is shown as overlay. 
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(c) Critical stage of topological transition 
with topological constraint. Diffuse wind- 
ing number is shown as overlay. 



(d) Past the critical stage of topological 
transition without topological constraint. 
Diffuse winding number is shown as over- 
lay. 
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(e) Past the critical stage of topological 
transition with topological constraint. Dif- 
fuse winding number is shown as overlay. 



Figure 3. Rotationally symmetric initial condition 



18 PATRICK W. DONDL, LUCA MUGNAI, AND MATTHIAS ROGER 

However, dE is not in the admissible class Ml as it cannot be parametrized by a 
single copy of 5*^. 

The reason why our topological constraint does not work properly in this example 
is due to the fact that Ti,{u) represents an approximation of the so called winding 
number, which depends on the orientation induced by E on the connected compo- 
nents of dE. In particular connected components of dE with opposite orientation 
(such as dBi/ii^) and 9Si/5((0, 3/4)) in the example above) compensate each other, 
and do not contribute to the value of the winding number (and consequently to the 
value of Ts{u)). A possibility to avoid such problem (firstly in the sharp-interface) is 
the following. 

Let (v?, r) denote a couple constituted by a finite collection F C -Bi(O) of W'^''^- 
regular, simple, closed and disjoint curves, and a function if G C^(T, [—1, 1]) such 
that iVrV'l ^ on r, that is (p assumes a constant value on each of the connected 
components of T. More precisely let A^ G N, 7^ G W'^^'^{S^, Bi{0)) {i = 1, . . . ,N) he 
diffeomorphisms such that (7j) fl (7^) = for i 7^ j, and let F = ((71), . . . , {'Jn)) and 
(f = d E [—1, 1] on (7j). We then set 

A{^,T):= f ^KrdV} = y^c, f k^JHK (A.2) 

Being 71, . . . 7Ar simple, regular, closed and disjoint curves, we can find (/i, . . . , In) G 
{1,2}^ such that, setting ip\r] := (—1)'' on (7j), we have 

t(F) := inf{A(<^,F) : ^ G BV{T), \VM{Q) = 0} = A{^[T],T) 

N 

Hence the functional T(F) counts the number of connected components of F, with- 
out taking into account of their orientation. It is then rather natural to look for a 
phase-fields approximation for T(F) in order to get a constraint on the topology of 
the diffuse interfaces stronger than the one obtained via Ti.{u). For this purpose we 
proceed as follows. We firstly consider a sequence A^{ip, u) of functionals defined on 
couples i^.u) G C^{VL, [—1, 1]) X C^(r2), and representing a diffuse interface approx- 
imation of y4((y9,F), Then, in analogy with T(F), we define a functional ^^(li) via 
minimization with respect to if of Afr{ip.,u). Finally we define the new topological 
constraint penalizing deviations of Tir{u) from 2tt. More precisely we start setting 



Ae{ip,u) := e^ / \V^\dxA / \Vu-^ ■V'^\e\Vu\dx 

-- [ (-eAu+-W'{u))^\Vu\dx. (A.3) 

A^{ip, F) formally presents a diffuse interfaces approximation of A{ip, F). The second 



term in ( A.3 ) represents a penalization of order e^'^ of the integral, with respect to the 
diffuse-length measure e\Vu\'^CiBi{o), of the variations of ip along the diffuse interface 
{y G -Bi(O) : Vu{y) 7^ 0}. Hence this term corresponds to a relaxation (at the diffuse 
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interface level) of the (sharp-interface) constraint iVrV^I = on F. The third term, 
as we have already seen in |4.22[ can be thought of as a phase- fields approximation of 
Jp (p K-p dh}. Finally the first term, whose contribution is infinitesimal of order e'^, is 
needed in order to ensure compactness in BV{VL, [—1, 1]) when minimizing As{ip,u) 
with respect to the first variable. In fact, we remark that, fixed u G C^(f2), we 
can apply the direct method of calculus of variations, and obtain the existence of a 
function Lp[u\ G BViVt, [—1, 1]) such that 

Ae{^[u],u)= inf A,((^,m), 
</3eci(Q,[-i,i]) 

where Ae{-,u) denotes the lower semi-continuous envelope of A^^-^u) with respect to 
the weak convergence in BV{VL^ [—1, 1]). Hence we define 



fe{u) 



J'riu) 



Co JBi(0) 



eAu + -W{u) ](p[u]\Vu\ dx, 



L] +e-^(fAu) 



27r 



(A.4) 
(A.5) 



and remark that when 99[m] = 1 the functional T^{u) coincides with the diffuse winding 
number TJu). 



In order to justify the choice of T^u) we first show (see Lemma A.l ) that if {u, 



is as in Lemma 
in Proposition 



4.1 



A.2 



%j£>0 



the value of T^iue) still converges to 27r as e — )• 0. Eventually, 
we analyze the behavior of Te(-) along sequences {us}e approx- 
imating (in an "optimal way") a finite collection F of simple closed, disjoint curves 
in -Bi(O), and obtain that T^Ue) converges to T(F). 



Lemma A.l. Let u^ he as in (4.10). Then we have 

Ym\fi,{ue) = 271. 

e-s>0 



(A.6) 



Proof. As in the proof of Lemma [41] we calculate that up to exponentially small term 

f (^-eAu, + -W'{u,))^\Vu,\ (A.7) 

Jb, (0) ^ 



eq[{tY K{s)Lp{'y{s) + tu{s)) dt ds. 



cos q;(s), sin a(s)) and observing that k{s) 



(A.8) 
a'{s) we deduce 



Writing 7'(s) 
that 

-L f-L 

K{s)^{-i{s)+tv{s))ds = I a{s)V^{-i{s)+tv{s))-{l + tK{s))-^'{s)ds 



We next observe that 



l\s) 



-<^(7(0)+tz/(0))(a(L)-a(0)). 



(A.9) 



Wu 



(7(s)+Ms)) 
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and define a{x) as the angle in [0, 27r) such that 



IVmI 



X 



cosal^a;^ 
sina(x] 



We then obtain for the first term on the right-hand side of (A.9) 

-L 



£q'^{tfa{s)V(p{-f{s) + tu{s)) ■ (1 + tK{s))y{s) dt ds 



a{x) 



{\d\<S/2} 






(x) ■ Vip{x)e\Vu\'^ dx\ 



<2n I \Vu^ ■ Vv^l^lVul 
'bi(o) 



(A. 10) 



For the second term on the right-hand side of (A.9) we have 



eq',{t)MliO) +HO)){a{L) - a{0)) 



< 27rco. 



(A.ii: 



This shows that for arbitrary ip G BV{Bi{0)] [—1, 1]), up to exponentially small terms 
in e 



Mv) > 



1 



Vm-^ • Vv? £:|Vm| -27r-27r 



Bi{0) 



Vm^ • Vv? £:|Vm| > -27r. 



Bi{0) 



On the other hand for 99 = 1 we have 

A,{y^) = -T,{u) ^ -277. 

This shows that lim^.^o^el'") = 27r, as claimed. 



D 



An application of the previous lemma shows that the improved topological con- 
straint in the case of a finite collection of simple, disjoint curves adds up the winding 
numbers of each single curve. In particular, these configurations are strongly penal- 
ized by the modified functional ^£. 

Proposition A. 2. Let E CC -Bi(O) be an open subset with C'^-bounday dE = 
^jLiilj) where '-fj are C^-diff'eomorphisms of the unit circle. Let {ugje C C^{Q) 



be constructed as in Lemma 4-i such that in particular Us — ?■ 2xe — 1 in L^(i?i(0)). 
Then 



limTpf-Up 



A^27r. 



e-s>0 



Proof. Since jumps of ip from 1 to —1 in the space between two connected components 



of dE are infitesimal of order e"', it is enough to repeat the proof of Lemma A.l for 
each connected component of dE. D 
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